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We investigate numerically the critical current of two-dimensional fully frustrated arrays of resis- 
tively shunted Josephson junctions at zero temperature. It is shown that a domino-type mechanism 
is responsible for the existence of a critical current lower than the one predicted from the transla- 
tionally invariant flux lattice. This domino mechanism is demonstrated for uniform-current injection 
as well as for various busbar conditions. It is also found that inhomogeneities close to the contacts 
makes it harder for the domino propagation to start, which increases the critical current towards 
the value based on the translational invariance. This domino-type vortex motion can be observed 
in experiments as voltage pulses propagating from the contacts through the array. 
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The two-dimensional (2D) Josephson junction arrays 
(JJA's) are ideal testing grounds for static and dynamic 
properties related to vortices rfnd constitute an intrigu- 
ing field of physics by itself.EJ They can be fabricated 
with high precision and are ideally suited for experi- 
ments and simulation studies, as well as theoretical ap- 
proaches!] Vortex physics is also of great importance for 
understanding the properties of high- T c superconductors, 
in particular, the interplay between vortex correlations, 
pinning, and dynamics. □ The present paper describes an 
interesting example of such an interplay. It focuses on 
what happens at the current threshold for the onset of 
resistance in the case of a 2D JJA in a perpendicular 
magnetic field. The threshold is in such a case caused 
by the onset of motion of the vortices induced by the ex- 
ternal magnetic field. We find that this onset of vortex 
motion is caused by a domino-type effect corresponding 
to voltage pulses traveling through the sample. 

The value of the critical currents in 2D JJA^s have 
been investigated theoretically by several authors U These 
studies have revealed that the zero-temperature critical 
current I c depends strongly on the value of the frustra- 
tion / defined by the number of flux quanta per plaquette 
induced by the external magnetic field, and they unan- 
imously have obtained the value I c = V2 - 1 « 0.41 
(in units of the critical current of a single junction) for 
the fully frustrated (/ = 1/2) square array. On the 
other hand, numerical calculations based on the resis- 
tively shunted junction (RSJ) model with the uniform- 
current injection method (see Ref. [| far-details) have per- 
sistently given the value I c = 0.35(1) Jj~l3 The discrepancy 
has usually been attributed to cither finite-size effects 
or an artifact of the uniform-current injection method, 
but without any really substantiating evidence for either 
claim. As to the first suggestion, the same numerical 
value has now been found in a very large array (256 x 256) 
(Ref. ||) and the absence of any size dependence up to this 
array size appears to preclude the possibility of an expla- 
nation in terms of finite-size effects. The second possibil- 



ity also by now seems less likely since attempts to model 
the injection of current through busbars, in close anal- 
ogy with actual experimental situations (see Ref. |Lfo« 
details), in fact seem to give an even smaller valueTffl'El 
Furthermore, in Ref. ^ an improved busbar method was 
proposed and I c » 0.35, close to the one obtained for the 
uniform-current injection, was found. Thus a value lower 
than the theoretically predicted one has been obtained 
with a variety of injection methods. 

In the present paper we resolve this long-standing dis- 
crepancy and find that the true critical current is in 
fact lower than the theoretically predicted one and cor- 
responds to a domino-type mechanism starting from the 
contacts which breaks the translational invariance of the 
vortex lattice. This means that the theoretically pre- 
dicted value, which presumes translational invariance, in 
fact gives an upper limit of I c . We conclude that the true 
critical current is always lower and that the precise value, 
because it is given by a domino mechanism, is sensitive 
to the inhomogeneities close to the contacts. 

We numerically study the fully frustrated L x x L y RSJ 
square arrays at zero temperature. We use the periodic 
boundary condition in the y direction and the external 
current Id is inserted on the left boundary (x = 0) and 
the extracted on the right boundary (x = L x ). The cur- 
rent Iij from site i to site j is the sum of the supercurrent 
and the normal current: 



Jij sm(9i 
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where J y - is the critical current of the single junction, the 
fully frustrated case corresponds to the magnetic bond 
angle A y - = ttx for all vertical links at x and zero for all 
links in the x direction, and finally, 6i — 8j is the normal 
current (in a convenient unit system). The RSJ equa- 
tions arej-pbtained by demanding current conservation on 
each siteo and are integrated either by starting from the 
ground-state configuration found from the simulated an- 
nealing Monte Carlo method for Id — or from a random 
configuration. The current injection methods are shown 
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in Fig. |j] and arenas follows: The uniform-current (UC) 
injection methodu corresponds to Jij — 1 for all links 
[Fig. |](a)]; the conventional busbar (CB) methods where 
the links on the boundaries have large corresponding 
to a superconducting busbar (we use Jy = 10 for the fai. 
links in Fig. [I]); Simkin's improved busbar (SB) methodtl 
where in addition to the superconducting busbar there 
are no magnetic fields on the plaquettes adjacent to the 
boundaries [Fig. |l|(c)]; the busbar method with inhomo- 
geneities near the boundaries (IB) [Figs. |l](d) and (e) are 
two variants, IB1 and IB2, respectively]. The idea be- 
hind these inhomogencities close to the contacts is that 
they create barriers for the onset of the vortex motion. 
The fat horizontal links in Figs. 1(d) and (e) play the 
role of such barriers. 

The translationally invariant ground state is a checker 
board pattern with a vortex at every second plaquette 
[see Fig. ||(a)]. Let us first assume the translational sym- 
metry even when an external current Id is turned on. 
Every second horizontal junction will then by symmetry 
have the same current and we below focus on the horizon- 
tal junctions with the largest current. For some external 
current Id — I s the current through these links reaches 
the maximum supercurrent for a single junction, which 
is 1 in our current unit. A straightforward calculation 
shows that this occurs at I s = (3 — \/5)/2 (see Table |). 
For Id < I s the ground state vortex configuration is al- 
ways stable. However, at some 1^ = J c > / s this ground 
state configuration becomes unstable and dissipation sets 
in. In the translationallyjiivariant case this happens at 
I c = \/2 - 1 (see Table 1 Jl2 which can readily be verified 
with the boundary condition in Ref. || which preserves 
the translational invariance; I c — 0.4142(1) has been ob- 
tained in excellent agreement with the prediction y/2— 1. 

In reality, however, the current pattern enforced by 
the contacts is not fully compatible with the translation- 
ally invariant vortex pattern. One consequence of this 
is that the external current I s b needed before a horizon- 
tal junction at the boundary reaches its maximum su- 
percurrent is lower than for the translationally invariant 
case. Table | gives the actual values for the boundary 
conditions shown in Fig. [j] together with the vortex col- 
umn for which it occurs. This lower value of I s b suggests 
that the vortex columns close to the contacts may start 
to move and dissipate energy before the rest of the vor- 
tices. Figure || illustrates some of the possibilities: The 
transition from Fig. ||(a) to ||(b) is the translationally 
invariant case where the whole vortex lattice moves in 
one piece. The transition from Fig. |(a) to |(c) followed 
by ||(c) to ||(a) illustrates the case when only the vortex 
column closest to the boundary moves and the others re- 
main fixed. Note that this is a boundary effect in the 
sense that it only produces a voltage drop for the hor- 
izontal junctions closest to the contacts and no voltage 
drop across the rest of the sample. Another possibility is 
that the motion of the first vortex column [Fig. ||(a) to 
H(c)] is followed by the second column [as illustrated by 



Fig. ||(c) to ||(d)], and then followed by the third and so 
on without stopping. Such a consecutive motion of vor- 
tex columns starting from the contacts and continuing 
without stopping we refer to as a domino-type motion. 
We have found that this domino-type motion starts at a 
lower current than what is required for the translation- 
ally invariant motion and that consequently the threshold 
current for the domino motion is the true critical current. 
Thus the domino effect is responsible for the critical cur- 
rent for all the boundary conditions shown in Fig. [j] and 
the corresponding values are given in Table |. 

The domino instability has two basic ingredients: The 
first is that the boundary and current injection break 
the translationally invariant current pattern close to the 
boundary which lowers the threshold current for the first 
vortex column to move. This is related to that I s \ for 
the first column is smaller than / s for the translationally 
invariant case (see Table |). The second ingredient is 
that this motion creates an instability which propagates 
through the whole sample. A hand-waving argument for 
this goes as follows: One first notes that as long as the 
checkerboard pattern of vortices is present a horizontal 
junction with lower current is always followed by one of 
higher current. However, when the first column moves 
[as illustrated by Fig. ||(a) to ||(c)] then two vortices 
become neighbors. This means that now two horizontal 
junctions with higher currents follow each other. Because 
of current conservation this means that the junctions in 
the second column with higher currents will momentarily 
further increase which makes this column unstable. The 
motion of the second column in the same fashion causes 
the third column to become unstable and so on. 

Figures |^(a) and (b) illustrate the domino mechanism 
for the uniform current injection: Figure ||(a) shows that 
the voltage is zero up to a critical current I c . The 
fact that the voltage per length saturates to a nonzero 
value above I c as the distance between the contacts is 
increased, demonstrates that it is a true onset of dissi- 
pation across the whole sample. Figure ||(b) shows the 
voltage drop V(x, t) at time t for the horizontal junc- 
tion at position x along a row of junctions in the pres- 
ence of a current slightly larger than the critical current 
[Id = 0.36 > I c = 0.350(1)]. As seen the voltage pulse 
always travels from one of the contacts towards the op- 
posite in accordance with the domino mechanism. In 
addition, the interference between domino propagations 
traveling in opposite directions creates interesting pat- 
tern in this 3D plot. A measurable consequence of this is 
that, if the voltage drop across different positions along 
the current direction is measured as a function of time, 
then traveling voltage pulses can be observed. This is in 
contrast to the translationally invariant case where the 
voltage appears at the same time along the whole current 
direction as is illustrated in Fig. ||(c). This translation- 
ally invariant case is simulated with the fluctuating-twist- 
boundary method, as is described in Ref. |[ 

The domino mechanism applies to all the situations 
shown in Fig. [j] with some variants: For example, in the 
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case of the conventional busbar injection [see Fig. |l|( 
the first vortex column becomes unstable at a current I c b 
but the domino avalanche does not set in until a larger 
current I c is applied (see Table Q). This is illustrated in 
Fig. [I], which shows the onset of dissipation at I s j, for a 
finite system and a further increase at I c . However, as 
the distance between the contacts becomes larger dissi- 
pation between I c b and I c decreases towards zero. The 
situation between I c b and I c is hence a boundary effect 
of the type illustrated by Fig. ||(a) to |2|(c) followed by 
||(c) to ||(a), i.e., only the first column moves. 

In Ref. [To] the critical current for a fully frustrated 
square array of Josephsoiv^unctions was measured and 
found to be I c w 0.42(2)Jb As apparent from Table | 
this value is larger than what is obtained for the RSJ 
model with uniform-current injection [Fig. |lj(a)] and bus- 
bar injections [Fig. [lj(b) an d ( c )] and is in fact close to 
the translationally invariant value. This means that the 
domino mechanism is somehow suppressed. One way of 
suppressing this mechanism is to introduce barriers close 
to the boundary which prevents the domino avalanche 
from spreading. The situation in Fig. |lj(d) shows an ex- 
treme variant of this and as seen from Table | this results 
in a critical current close to the translationally invariant 
case. However, also in this case the domino mechanism is 
responsible for the onset of dissipation. Thus the transla- 
tionally invariant value appears to be an upper limit and 
the true critical current is always lower and given by the 
domino mechanism. A perhaps more realistic situation 
is with some imperfections close to the boundaries which 
creates barriers for the vortex motion like in Fig. ^(e) and 
as seen from Table | the critical current is again close to 
the translationally invariant value. Figure [sj shows the 
domino mechanism for this case in a 3D plot. [Note the 
similarity with Fig. |^(b).] 

The time scale of the Josephson junction array is given 
by to = TijleR^J with the shunt resistance Rn- For the 
1000 x 1000 SNS junction array in Ref. M where R N m 2 
mfi and J w 7 mA this means to » 10"^ sec. The time 
for a voltage pulse to reach the middle would then for 
the uniform current injection [Fig. 3(b)] as well as for 
the case with inhomogeneities (Fig be of the order 
of 10 -7 sec. Thus it seems that it should be possible 
to verify the domino mechanism by either direct voltage 
measurements or by some vortex imaging technique.E3 

In summary, on the basis of simulations for the RSJ 
model, we have concluded that the critical current for a 
2D fully frustrated Josephson junction array is due to a 
domino mechanism and that the voltage pulses resulting 
from the domino avalanches should be open to experi- 
mental verification. 
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TABLE I. The critical current (7 C ), the boundary critical 
current for the motion of the first vortex column (Icb), and 
the external current needed for the first junction to reach its 
maximum current (I s b) (in hi, i— 1,2,5 denotes the column it 
occurs in). TI denotes the translationally invariant case; see 
Fig. [j] for the other current injection methods. 

7c Icb Isb 

TI y/2 - 1 « 0.4142 I ab = Is = (3 - V5)/2 w 0.382 

UC 0.350(1) Li = 0.285(1) 

CB 0.324(1) 0.277(1) I s i = 0.183(1) 

SB 0.337(1) L2 = 0.265(1) 

IB1 0.4135(5) Iss = 0.366(1) 

IB2 0.4138(5) 0.281(1) I s i = 0.249(1) 
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(a)UC 

□ f=m 

□ f=o 
I 10/ 



(b) CB (c) SB 
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FIG. 1. Current injection methods: (a) Uniform current 
(UC) injection method, (b) conventional busbar (CB), (c) 
busbar method in Ref. W where unfrustrated columns (/ = 0) 
are used near busbars (SB), (d) and (e) show busbar methods 
with inhomogeneities near boundaries (IB1 and IB2, respec- 
tively). The fat links have a much larger coupling strength 
than the normal links. 
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FIG. 2. The configuration with translational symmetry in 
(a) can evolve in different ways: From (a) to (b), the transla- 
tional symmetry is preserved by the rigid motion of vortices. 
From (a) to (c), only one column near boundary moves. Once 
the motion (a) to (c) happens, this unstable vortex configu- 
ration can be resolved by the motion either to (d) or to (a). 
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FIG. 3. (a) Time averaged voltage V per junction versus 
external current Id for uniform-current (UC) injection case. 
The absence of finite-size effects for large sizes shows that 
I c = 0.350(1) is the true critical current, (b) Voltage V(x,t) 
at time t across the horizontal junction at position x along a 
row of junctions for UC at Id = 0.36: voltage pulses travel 
from the contacts into the sample, (c) V(x,t) for the trans- 
lationally invariant case at Id = 0.417: the voltage occurs in- 
stantaneously across the whole sample corresponding to the 
rigid motion. 
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FIG. 4. Time averaged voltage V per junction versus exter- 
nal current Id for for the conventional busbar case [Fig. 0(b)] : 
V between the boundary critical current I s t ~ 0.28 and the 
true critical current 7 C « 0.32 vanishes as the distance be- 
tween the contacts increases because it is due to the motion 
of the vortex columns closest to the contact. The onset of 
nonzero voltage across the sample requires the domino mech- 
anism which sets in at I c « 0.32. 




FIG. 5. Voltage V(x,t) at time t across the horizontal 
junction at position x along a row of junctions for the IB2 
case with inhomogeneities close to the contacts in Fig. ^(e) at 
Id ~ 0.416. The critical current is in this case close to the up- 
per limit given by the translationally invariant case, but the 
domino mechanism is still responsible for the onset of nonzero 
voltage. 
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